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Abstract 

This note presents a self-contained and streamlined exposition of chemical-network results due to 
Feinberg, Horn, and Jackson. As an application, it shows the global asymptotic stability to equilibria 
of McKeithan's kinetic proofreading model for T-cell receptor signal transduction. 

1 Introduction 

This work was motivated a question which arose during Carla Wofsy's series of talks [Q. Consider the 
following system of first-order ordinary differential equations, for nonnegative functions Ci{t): 

Co = ki (T* - £f =o a) (M* - £f =o a) - (fc-i,o + fcp.o) Co 
Ci = kp i^ iCj_i — (fc_ i,i + Gi 

C N — kp N ^\C N —i k—i N C N 

(dots indicate derivatives with respect to time t) where the subscripted fc's, as well as M* and T*, are 
arbitrary positive constants. 

These equations represent the dynamics of the "kinetic proofreading" model proposed by McKeithan 
in H in order to describe how a chain of modifications of the T-cell receptor complex, via tyrosine phos- 
phorylation and other reactions, may give rise to both increased sensitivity and selectivity of response. 
The quantities C(i) represent concentrations of various intermediate complexes, and the assumption is 
that recognition signals are determined by the concentrations of the final complex C N . The constant k\ 
is the association rate constant for the reaction which produces an initial ligand-receptor complex Co 
from a T-cell receptor (TCR) and a peptide- major histocompatibility complex (MHC). The constants 
k Pi i are the rate constants for each of the steps of phosphorylation or other intermediate modifications, 
and the constants fc-i^ are dissociation rates. The differential equation for Co could also be written in 
an alternative manner, as 

C = fciTM - (A_i,„ + kp,a) C , 

where T(t) and M(t) represent the concentrations of TCR's and MHCs respectively. The two conser- 
vation laws 

M + Y,Z Ci = M*, T + ^f =0 Ci = T*, 
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when solved for T and M, and substituted back into the equation, give the form shown earlier. 

McKeithan's paper focused on the analysis of equilibria of the above differential equations (and made, 
for simplicity, the assumption that k Pt i = k p and fc— i,i = for some fixed k p and fc_i). The question 
that arose during |I| was: what can be said about the dynamics of these equations? 

We show here that the best possible conclusion is true: there is a unique equilibrium point, and this 
equilibrium is globally asymptotically stable. 

It turns out that this result can be derived as an almost immediate corollary of the beautiful and 
powerful theory of deficiency zero chemical reaction networks with mass- action kinetics developed by 
Feinberg, Horn, and Jackson, cf. jj], ||. In this note, we wish to provide a totally self-contained and 
streamlined presentation of the main theorems for such networks. In order to keep matters as simple 
as possible, however, we specialize to a subset ("single linkage class networks") which already contains 
the kinetic proofreading model; the reader is referred to 0] , and the bibliography therein, for analogous 
results concerning more arbitrary deficiency zero networks. 

The organization of this note is as follows. The results in question can be explained in terms of 
a special yet very general and appealing class of dynamical systems, which we introduce in Section |^, 
where the main theorems are also stated. In Section [| we study equilibria and invariance properties for 
this class of systems. Section [| has proofs of the stability theorems. 

In Section^, we specialize to the kinetic proofreading example. We do not actually use any terminol- 
ogy from chemical network theory, nor do we define the terms "deficiency zero" and "chemical reaction 
network," but also found in Section || are brief remarks concerning motivations for the general form of 
the systems considered, and the relation to chemical network concepts. 

We wish to emphasize that this note is mainly expository, and credit for the mathematical results 
lies with Feinberg, Horn, and Jackson. It is our hope that this exposition will serve to make a wider 
audience in the dynamical systems and control theory communities aware of their work. A follow-up 
note Q will add new results concerning robustness, dependence on parameters, and control-theoretic 
properties of the systems studied here. 
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2 Definitions and Statements of Main Results 

Some standard notations to be used are: 

• K>o (resp., R+) = nonnegative (resp., positive) real numbers 

• (resp., R™ xm ) = n-column vectors (resp., m x m matrices) with entries on R+; similarly for 

• Kg — boundary of R> , set of vectors x £ R> such that Xi = for at least one i £ {1, . . . , n} 

• x 1 = transpose of vector or matrix x 

• \x\ = Euclidean norm of vector in R™ 

• (x, z) = x' z, inner product of two vectors 

• V 1 - = {x \ (x,z) = OVz £ D}. 
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The systems to be studied are parametrized by two matrices A and B with nonnegative entries, as 
well as a nonnegative function 6, and have the following general form: 

m rn 

x = f(x) = ^^a ij e(xi)^e(x^...9(x n )^(b i -b j ) (1) 
»=i j=i 

where &^ denotes the £-th column of B (notice that the diagonal entries of A are irrelevant, since 
bi — bi = 0). Several restrictions on A, B, and 6 are imposed below. The powers are interpreted as 
follows, for any r, c > 0: r° = 1, C = if c > 0, and r c = e clnr if r > and c > 0. 

The main example of interest is when 6(y) = |y| and £? is a matrix whose entries are nonnegative 
integers. In that case, the equations (|IJ) are polynomial for nonnegative vectors x. However, we will 
allow a more general class of systems. 

We next describe the hypotheses on 9, A, and B. The map 

9 : K — > [0, oo) 

is locally Lipschitz, has 0(0) = 0, satisfies J |ln0(y)| dy < oo, and its restriction to R>o is strictly 
increasing and onto. We suppose that 

A = ( aij ) £ R^ l xm is irreducible (2) 

(that is, (/ 4- A)" 1 ^ 1 £ R™ xm or, equivalently, the incidence graph G(A) is strongly connected, where 
G(A) is the graph whose nodes are the integers {1, . . . , m) and for which there is an edge j — > i, i j, 
if and only if a,j > 0), and that 

B = (6i, . . . , 6 m ) e M>o m has rank TO ( 3 ) 
(so, its columns 6j are linearly independent), 

no row of £? vanishes, (4) 

and 

each entry of B is either or > 1 . (5) 

This last hypothesis insures that f(x) in (Q) is a locally Lipschitz vector field. 

From now on, we assume that all systems ^) considered satisfy the above assumptions. 

Our study will focus on those solutions of ([!]) which evolve in the nonnegative orthant K> . Recall 
that a subset S C R ra is said to be forward invariant with respect to the differential equation x = f(x) 
provided that each solution x(-) with x(0) S S has the property that x(t) G S for all positive t in the 
domain of definition of x(-). We show in Section || that the nonnegative and positive orthants are forward 
invariant: 

Lemma 2.1 Both R" and R™ are forward-invariant sets with respect to the system (Q). 

We will also show in Section || that there are no finite explosion times: 

Lemma 2.2 For each £ 6 R> there is a (unique) solution x(-) of (0) with x(0) = £, defined for all 
t > 0. 



In order to state concisely the main results for systems ([]]), we need to introduce a few additional 
objects. The subspace 

D := span {bi - bj, i ^ j} = span {bi - 6 2 , . . . , h - b m } (6) 
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can be seen as a distribution in the tangent space of K"; it has dimension to — 1 because adding b\ to 
the last-shown generating set gives the column space of the rank-m matrix B. For each vector p £ R", 
we may also consider the parallel translate of D that passes through p, i.e. p + T> ~ {p + d, dG ©}. A 
set S which arises as an intersection of such an affine subspace with the nonnegative orthant: 

S = (p + D)f|R5 

(for some p, without loss of generality in K" ) will be referred to as a class. If S intersects the positive 
orthant R™ , we say that S is a positive class. The significance of classes is given by the fact that any 
solution x(-) of ([!]) must satisfy 



5 (t) 



„t mm 

o(0) = / f(x(s))ds = ^^7(*)ft-y 62 ; 



where 7(4) = f* a ij 9(x 1 (s)) bl i 9(x 2 (s)) b ^ . . . 9(x n (s)) K * ds, so x(t) e au(0) + CD for all t. In particular: 

Lemma 2.3 Each class is forward invariant. □ 

We denote by E (respectively, E + or E ) the set of nonnegative (respectively, positive or boundary) 
equilibria of ([!]), i.e. the set of states x € K™ (respectively, £ K™ or 6 Kg), such that f(x) = 0. Of 
course, E is the disjoint union of E + and Eq. 

Theorem 1 Consider any system under the stated assumptions. For every maximal solution of (J^j 
with x(0) G R> ; ^ ftoZrfs </ia< a;(i) — > E as t — > +00. 

This will be proved in Section The invariance of classes (which are contained in subspaces of 
dimension m — 1 < n) precludes asymptotic stability of equilibria. The appropriate concept is that of 
asymptotic stability relative to a class. We say that an equilibrium x 6 S is asymptotically stable relative 
to a class S if it is (a) stable relative to S (for each e > 0, there is some S > such that, for all solutions 
x(-), \x(0) — x\ < S and x(0) £ S imply \x(t) — x\ < e for all t > 0) and (b) locally attractive relative 
to S (for some e > 0, if |x(0) — S| < e and x(0) € then x(t) — > x as f — > +00). We say that x E S 
is globally asymptotically stable relative to a class S if it is stable relative to S and globally attractive 
relative to S (x(t) — > 5 for all solutions with x(0) S 5 1 ). The main results are as follows; the first part is 
shown in Section ^ and the remaining two in Section ^. 

Theorem 2 Consider any system (jj[), under the stated assumptions. Fix any positive class S. 

a. There is a unique equilibrium xs 6 Sf]E+. 

b. The equilibrium x$ is asymptotically stable relative to S . 

c. The equilibrium xs is globally asymptotically stable relative to S if and only if Sf]Eo = 0. 



Example 2.4 The following trivial example may help in understanding the above theorems. We take 
n = m = 2, 9(y) = \y\, A = identity matrix, and 




The system (|I]) is (for nonnegative states): 

X\ = (x\ - 1)X\X2 
x 2 = 



4 





X 


»S ■■)■■ 
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Figure 1: Example 2.4, dark lines indicate equilibria 



and thus Eq = Rq = {x \ x\X2 = 0} and E + = {x \ x\ = 1,X2 > 0}. The positive classes are the sets 
S = S r = {x | x\ > 0, X2 — r}, for each r > 0, and for each such S = S r , x — (1, r)' is asymptotically 
stable with domain of attraction {x \ x\ > 0, X2 = r}. See Figure |]. Each class S r has a second 
equilibrium (0, r)', but this second equilibrium is in the boundary, so there is no contradiction with part 
a of Theorem |^. Regarding Theorem [l], observe that every trajectory either converges to an interior 
equilibrium (1, r)' or it is itself a trajectory consisting of an equilibrium (and hence also converges to E, 
in a trivial sense). □ 



2.1 Other Expressions for the System Equations 

The equations ([!]) have a considerable amount of structure, and various useful properties are reflected 
in alternative expressions for the system equations. 

Let us introduce the map 

p(y) := ln%) 

(with p(0) = —oo). The restriction 

p:R + ^ M 

is locally Lipschitz, strictly increasing, and onto R. It also satisfies lim^o p{y) — —oo and \p(y)\ dy < 
oo. For any positive integer n, we let 

p : W 1 -» [-oo, oo)™ : x i-> (p(x x ), . . . , p{x n ))' 

(we do not write "p n " to emphasize the dependence on n, because n will be clear from the context). 
Then ([!]) can also be written as 

m m 

x = f{x) = Y,Y, a n e{bi,m) ^ - h o) ■ ( ? ) 

i=i j=i 

Here, the expression "e^'^ 2 ^" in (Q) is interpreted in accordance with the conventions made for powers: 
if X is a vector and k £ {1, . . . , n} is an index such that Xk = and bkj > 0, then e bk ^ p ^ Xk ^ — 0, consistently 
with e~°° = 0, and thus also 

e (bj<P( x )) = e bljp{xi) e b 2 jp(x 2 ) e bnjp(x n ) _ q 

but, if = 0, then we have e bkjP ^ Xk ^ = 1. 

Another useful way of rewriting (]l|) is as follows. We write for the fc-th coordinate of / (i.e., the 
coordinates Xk of solutions x satisfy ik — fk(x)). The terms in the sums defining fk can be collected 
into two disjoint sets: those that do not involve a product containing 9(xk), for which bkj = 0, and 
those which do involve 9(xk)- The latter, by assumption (|b|), have bkj > 1, so we can factor 0{xk) 
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from 6(xi) hl: >6(x2) b2j ■ ■ ■ 8(x n ) bn ^ and there remains a locally Lipschitz product. In other words, we can 
introduce, for each k G {1, . . . , n}, these two locally Lipschitz functions: 

a k (x) := (f>« {b M -b kj )\e(x 1 ) b »0{x 2 ) b ^...0{x k ) b »- 1 ...0{x n )^ (8) 

jeJ fc ,i \i=i / 

and 

Pk(x) := Yl (E«u^)^i) 6lJ ^2) 62j '...^„) 6 -, (9) 
je.7 fe ,o \»=i / 

where J k ,i ■= {j \ b k j > 1} and J kfi ■= {] I b k j = 0}. In terms of these, 

fk(x) = a k {x)6{x k ) + p k {x). (10) 

In particular, since 0(0) = 0, 

*k = => f k (x)=p k {x) (11) 

so, since /3k(x) > for all x, the vector field / always points towards the nonnegative orthant, on the 
boundary Rg . 



3 Equilibria and Invariance 

We start with a global transversality result for D and D^. Then, we study interior and boundary 
equilibria, and provide basic results concerning the behavior of nonnegative solutions. 

3.1 A Coordinatization Property 

The next result shows provides, when specialized to classes, a one-to-one correspondence between points 
in R™ and pairs (<S, 1Z) consisting of a positive classes S = p + S and "coclasses" 1Z = p~ 1 (q + p(S ± j). 

Lemma 3.1 Let D be any subspace of W\ For each p, q in R™ , there exists a unique x = tp(p, q) £ R™ 
such that: 

x-p e D (12) 

and 

p{x) - p{q) e . (13) 
Proof. We start by introducing the following mapping, for each i £ {1, . . . , n}: 

rt+p(qi) 

Li(t) := / P^ 1 (s) ds - pit 

Jo 

defined for ( e I. Since p _1 is an increasing onto map from R into R+, L'^t) = p _1 (t + p{qi)) — Pi, 
and hence also Li(t), increases to infinity as t — > +oo. Also, Li(t) — > +oo as t — > — oo, because is 
nonnegative and pi > 0. Thus, Li is proper, that is, {t \ Li(t) < v} is compact for each v. 

Now we take the (continuously differentiable) function 

n 

Q(y) := ^LifoO 

i=l 

thought of as a function of y G R™. This function is also proper, because 

n 

{y I Q(y) < w} c ]J{t | Li(t) <w-(n- 1)1} 

i=l 



G 



where £ is any common lower bound for the functions Li. Restricted to D, Q is still proper, so it attains 
a minimum at some point y E D. In particular, y must be a critical point of Q restricted to D, so 

(VQ(y))' = {p- 1 {y 1 + p(q 1 ))-p 1 ,...,p- 1 (y n + p(q n ))-Pn)' G D ± . (14) 
Pick x E R+ such that p(x) — y + p(q). Then p(x) — p(q) E D by definition, and ( |T^ ) gives also 

Finally, we show uniqueness. Suppose that there a second z E K™ so that z~p £ D and p(z) — £ 
D . This implies that x — z £ D and p(x) — £ £> . Since p is an increasing function, we have that, 
for any two distinct numbers a, b, (a — b)(p(a) — p{b)) > 0. So 

n 

}^Xj - Zi)(p(xj) - p{zj)) = (x - z,p(x) - p{z)) = 0. 

i=l 

implies x — z. I 
The following quantity measures deviations relative to D . Let us define, for each x, z € R™ : 

mm 2 

Note that 5 (a;, z) = if and only if p(x) — p(z) £ D^, since D is spanned by the differences &j — 6j, 

Remark 3.2 We could also have defined a smaller, but basically equivalent, sum using only the gener- 
ating differences &j — b\, i = 1, . . . , m — 1, but the above definition for 6 seems more natural. Moreover, 
note that if we let 

m— 1 n 

A(x, z) := i( b i ~ b i,p{x) - p(z))f + Y (( v ti x ~ z )f > 

i— 1 i—m 

where the i>j constitute a basis of D^, then the uniqueness part of Lemma 3.1, applied with D = D, 
gives that A(x, z) — if and only if x = z. (Because x = <p(x,x), and A(x, z) = implies z = ip{x,x).) 
□ 



3.2 Equilibria 

It is convenient to also express the dynamics (mj in matrix terms. Letting 



A = A- diag ^ a iX , . . . , ^ < 



021 - a ^2 



V 



a>n2 



we write: 



where b is the mapping 



x = f{x) = BA 8s(x) 



pn . Tram 



p(*)> 



air. 

«2, 



Em . 



(16) 



obtained as the composition of the maps x h- ► p(x), z i— > B'z, and y h- ► (e Vl , . . . ,e Vm )' . In particular, 
since rank B = m and p maps M_|_ onto R, 



the restriction Qb 



Dn , irpm 



is onto. 



(17) 
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Note that e B (a;) = {0(x 1 )^6(x 2 ) bai ■ ■ ■ 9{x n ) b ^,. . . , 9{x 1 ) b ^8{x 2 ) b ^ . . . 0{x n ) b ™ )'. 
For any two x, x € R™ , it holds that: 

(3 k > 0) Q B (x) = kOb(x) <=^> p(x) - p(x) e D 1 - (18) 

(recall the definition (^) of D). To see this, denote y :— Qb(x), y :— 8b(x). If j/ = ny, then, with 
fc := ln/t, lnj/j = + lnj/j, for each j = 1, ... ,m. Thus (bj,p(x)) = k + (bj,p(x)), which implies that 
(bj,p(x) — p(x)) — k for all j, and therefore 

(bj - bi, p(x) - p(x)) = Vi,j. 

Conversely, if this holds, we may define k := (bi, p(x) — p(x)), K := e k , and reverse all implications. 

We also note, using once again that B has full column rank, that f(x) = is equivalent to AQb(%) = 
0, that is: 

Lemma 3.3 A state x is an equilibrium if and only if Qb{x) G kcr A. □ 

We now consider the matrix A. The row vector 1= (1, . . . , 1) has the property that 1A = (0, . . . , 0), 
so, in particular, A is singular. The following is a routine consequence of the Perron-Frobenius (or finite 
dimensional Krein-Rutman) Theorem. 

Lemma 3.4 There exists y <E R™ |~| ker A so that (R£ \ {0}) f| ker A = {ny, k > 0}. 

Proof. If y G R.> is any eigenvector of A, corresponding to an eigenvalue A, it follows that = lAy = 
l\y = Xq, where q := ly is a positive number (because y, being an eigenvector, is nonzero), and therefore 
necessarily A = 0. In other words, a nonnegative eigenvector can only be associated to the zero eigenvalue. 
Pick now any 7 > large enough such that all entries of A := A+^/I are nonnegative. Since the incidence 
graph G(A) coincides with G(A), it follows that A is also irreducible. By the Perron-Frobenius Theorem, 
the spectral radius a of A is positive and it is an eigenvalue of A of algebraic multiplicity one, with an 
associated positive eigenvector y G R" . Moreover, every nonnegative eigenvector y G R™ associated to 
a is a positive multiple of y. As adding 7/ moves eigenvalues by 7 while preserving eigenvectors (that 
is, (A + "fl)y — (A + 7)2/ is the same as Ay = Xy), y is a positive eigenvector of the original matrix A. 
It is necessarily in the kernel of A, since we already remarked that any nonnegative eigenvector must be 
associated to zero. Finally, if y is any other nonnegative eigenvector of A, and in particular any element 
of (R> \ {0}) Piker .A, then it is also a nonnegative eigenvector of A, and thus it must be a positive 
multiple of y, completing the proof. I 



Corollary 3.5 The set of positive equilibria E+ is nonempty. Moreover, pick any fixed x G Ej< 
for any positive vector x G R™ , the following equivalence holds: 



x G -Ea 



5(x, z) = . 



Then, 



(19) 



Proof. By Lemma 3.4, there is some y G R™ in ker A. By (|l7|), there is some x G R™ such that &b{x) = y. 
In view of Lemma 3.3, x is an equilibrium. 

Now fix any x G E + and any x G R™, and let y := Qb(x), y := Qb(x). Suppose x G E + . By 
Lemma 3.3 , y G ker A. By Lemma [3.4| , every two positive eigenvectors of A are multiples of each other, 
so there is some k G R+ such that y — ny. By (|l8|), p(x) — p(x) G D^. Conversely, if this holds, then, 
again by (|lq), y = ny. hence y is also an eigenvector of A, so by Lemma 3.3 we conclude x G E + . I 
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3.3 Proof of Part a in Theorem [2] 



Pick any positive class S = p+D, p G R™ . By Corollar y |3.5| , there is some equilibrium x G E + . We apply 
Lemma [O] with D = D, to obtain xs = (f(p,x). By (|13D and (|l9|), xs £ £'+. By (|l^), - p G B, i.e., 



also xs G S, as required. To show uniqueness, suppose that also z £ 5*p| Since z £ S 1 , — z G D, 
and since z £ £+, pl^s) — p(z) G D 1- (by Corollary |l^) . Thus the uniqueness assertion in Lemma 3.1 
gives xs = z. I 

3.4 Boundary Equilibria 

Fix any x = (x\, . . . , x n ) G K> - We wish to study the implications of some coordinate Xk vanishing. 
For each j G {1, . . . , to} we consider the set 

Sj : {/.' | b kj > 0} 
which is nonempty, by (||). Note that fc G Sj if and only if j G Jk,i- 

We will use repeatedly the following fact, for any j G {1, . . . , to} and fc G {1, . . . , n}: 

x k = and k G Sj => 6l(a;i) bl ^(x2) b2j . . . <9(x„) b " J = (20) 
which is obvious, since 9(xk) bkj vanishes when Xk = and bkj ^ 0. In particular, 

x k = 0(x 1 ) bl ^(a; 2 )^ . . .0(a;„) h -6 fe , = (21) 
since either 6/-j = or (|2^) applies. 

Lemma 3.6 Take any i, j G {1, . . . , m} such that a.y ^ 0. Then, 

(V^ G Sj) [xt > 0) => (Vfc G §i) (x fc > or f k (x) > 0) . (22) 

Proof. Pick any fc G Si, and assume that xj, = 0. The assumption "£ G Sj => xg > 0" is equivalent to 
j G J^i > 0. So, since Xfc = 0, necessarily j G Jfc,o- By (|TT|) , fk{x) — (3k(x), where (3 is as in (||), 

and the index j being considered does appear in the sum defining Moreover, for each £ G Sj, xt ^ 
by hypothesis, so ^xi) 6 ^^) 6 ^ . . . 9{x n ) bnj > 0. On the other hand, k G Si means that bki > 0, and 
also 7^ 0. Thus, the term involving this particular i and j in the sums defining [3k(x) is positive (and 
the remaining terms are nonnegative) . I 

Proposition 3.7 Take any x G IR™ . The following properties are equivalent: 

1. x G £o- 

2. For every j G {1, . . . , to} there is some k G Sj such that Xk = 0. 

3. For every j G {1,...,to}, 6>(a;i) fc « (9(a:2) b2j ...9{x n ) b ^ =0. 

Proof. [1 2] Pick any a; G E . If the second property is false, then there is some index j such that 
xg > for all £ E Sj. We claim that for every index i, Xk > for all fc G Si. Since (J . Sj = {1, . . . , n} 
(recall hypothesis (Q)), this will mean that Xfc > for all fc, so x could not have been a boundary point, 
a contradiction. Let J = {j \ xg > V£ G Sj} and let / = {1, . . . , m} \ J. We know that J ^ and must 
prove that 7 = 0. Suppose by contradiction that J ^ 0. Pick some i £ I and j & J such that aij 7^ 



(irreducibility of A), and take any fc G Si. Lemma 3J: gives that either x k > or fk{x) > 0. Since a; is 



an equilibrium, fk(x) = 0. So Xfe > for all fc G Si, contradicting the fact that i G I. 

[2 => 3] Suppose that the second property holds. Pick any j G {1, . . . , m}. As there is some fc G Sj 
such that x k = 0, © gives 6{x x ) b v 9{x 2 ) b ^ . ■ ■ 9{x n ) b ^ =0. 

[3 =S> 1] Obvious. I 
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3.5 Invariance Proofs 



The key technical fact is as follows. 

Lemma 3.8 Suppose that x : [0,£*] — > R" is any solution of ([!]) such that x(t) G R> for all £ G [0,£*]. 
Then the following implication holds for any k — 1, . . . , n: 

x fe (o)>o => x fe (r)>o. 

Proof. Suppose that k is so that Xk(0) > 0. We consider, for £ G [0,£*], the scalar function y(t) := X}~(t) 
and the functions a(t) := otk(x(t)) and (3(t) := fik(x(t)). By (fio|), 

y(t) = «(£)%(£)) +/?(£). 

Since is locally Lipschitz, j3(t) > for all £, and y(0) > 0, it follows by a routine argument on differential 
equations that y(t*) > 0, as desired. (The argument is: Let z solve z — a(t)9(z), z(0) = y(0) > 0. Since 
9 is locally Lipschitz, and is an equilibrium of this equation, z{t) > for all £ in its domain of definition. 
Moreover, we have that z = gi(t, z) and y = gi(t, y) with gi(t,p) < g2(t,p) for all p (because f3 > 0). By 
a standard comparison theorem, we know that z(t) < y(t) for all £ in the common domain of definition 
of z and y. Since y(t*) is well-defined, z(t) remains bounded, and thus is defined as well for t = t* . So, 

y(t*) > z(t*) > o. ■ 



Remark 3.9 Note that the only facts used in the proof are that 9 is locally Lipschitz and that 9(y) > 
for y > (which implies fi > 0); 9(y) > for all y is never needed. We assumed 9 > just for convenience 
when displaying the general form of the systems being considered. □ 



Corollary 3.10 The set R™ is forward invariant for ([!]). 

Proof. Consider any solution x : [0,T] — > R™ of (|l|), and suppose that x(0) € R™. We must prove 
that x(T) £ R™. Since x(0) is in the interior of R> , the only way that the conclusion could fail 
is if x{t) £ Rq for some t G (0, T]. We assume that this happens and derive a contradiction. Let 
t* := min{i G [0,T] \ x(t) G Rq} > 0. By minimality, for all i, x t (t) > for all t G [0,t*), and in 
particular Xi(t) G R" for all £ G [0, £*], and also there is some index k such that Xk{t*) = 0. But this 
contradicts the conclusion of Lemma [T^. I 

The closure of an invariant set is also invariant, so: 
Corollary 3.11 The set R™ is forward invariant for (nl). 



Lemma 3.12 Consider any solution x : [0, T] — > R™ of (|l|) for which x(0) G R> - Suppose that there is 
some j G {!,..., to} such that x e {0) > for all I G §j. Then, ac(t) G R^ for all~£ G (0, T\. 



Proof. We know that x(t) G R> , by Corollary 3.11. We need to see that every i belongs to the set 



I := {i G {l,...,TO}|a; fc (£) > OVfc G S,,V£ G (0,T]} 
since then lj i §i — {1, . . . , n} gives the desired conclusion. Pick any i such that ^ and any k G 



Then, Lemma says that Xk(0) > or A-(x(0)) > 0. If Xk(0) > 0, Lemma |3.8| , applied on any 
subinterval [0,£*], says that Xk{t) > for all t. If, instead, fk{x(0)) > 0, then ifc(0) > and Xk(0) > 
imply that Xk{t) > for all £ small enough, so also (again by Lemma |3~8| ) for all £. Thus i E I, which is 
therefore nonempty. 
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Suppose that H := {1, . . . , m} \ J ^ 0. Let s € J and /i £ ff be so that a/ji ^ (irreducibility of 
A). We will show that, for any given to G (0, T], and for any given k G S/j, Xk(to) > 0, and this will 
contradict h G P . 



3.£, which now says 



3.8 that x fe (£o) > 0. 



Since i G J, Xi(to/2) > for all £ G S,;. Then, we can apply once again Lemma 
that Xk{t /2) > or f k (x(t /2)) > 0. As before, Xk(t /2) > implies via Lemma 
And if /fc(x(t /2)) > 0, then i fe (t /2) > and Xk(to/2) > imply again that x k (t a ) > 0. I 

Corollary 3.13 Consider any solution x : [0, T] -> K£ of (g) for which x(0) £ P - Then, x(t) G M™ 
for allf G (0,T]. 



Proof. By Proposition 3.7, x(0) G' Po implies that there is some j G {1, . . . , to} such that ^ for all 



^ G Sj. So, Lemma 3.12 insures that x{t) G K™ for all t G (0, T 



We conclude that every trajectory starting on the boundary Kg which is not an equilibrium must 
immediately enter the positive orthant. 

4 Stability Proofs 

We start by establishing some useful estimates. 
Lemma 4.1 Define the following quadratic function: 

rn m 

Q(Vl, ■■■,Vm) ■= X] X] ai ^ 7n ~ ^ • ^ 23 ) 
i=l j=l 

Then, there exists a constant k > such that 

m m 

Q(qi,...,q m ) > ^^(«-9 3 ) 2 (24) 
i=l j=l 

for all (q u ...,q m ) eW m . 

Proof. We first observe that 

Q(qi, ...,q m )=Q => qi = q m , i = 1, • ■ ■ , m - 1 . 

Indeed, obviously Q(qi, ■ ■ ■ , 9m) = implies q t — qj for each pair i,j for which et.y ^ 0. Now let / be the 
set of indices i such that qi = q m , and J its complement; as m G J, J 7^ 0. We need to see that J = 0. 
Suppose that J 7^ 0. The connectedness of the incidence graph of A provides an i G / and j G J such 
that ay 7^ 0. Thus, q.j = q% = q m , contradicting j G J. 

Now consider the following quadratic form in m — 1 variables: 

m — 1 m — 1 m — 1 m — 1 

P(£i, . . . ,Cm-i) - X! E fl y ~ ^') 2 + X flim ^ 2 + X a ™i£j ■ 

Since - ?7 m ) - (77.,- - 7j m ) = 77,; - 77,- for all z, j, one has 

Q(r)l,---,Vm) = P(Vl-Vm,---,r)m-l-Vm)- 

Note that P is positive definite: if P(qi, . . . , <7 m -i) = 0, then Q(qi, . . . , q m -i, 0) = 0, which as already 
observed implies that all qi = 0. Thus, there is some constant kq > such that 

m — 1 

2 



P(pi, . . . ,/Jm-l) > K X P ' 
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for all (pi, . . . ,p m -i) £ K. m 1 , which means that 

m— 1 

Q(<?i, ■ ■ ■ , q m ) > k E - <? m ) 2 (25) 



j=i 



for all (qi, . . . , q rn ) G E m . As (qi — qj) 2 < 2(q, — q m ) 2 + 2(qj — q m ) 2 for all i, j, we may re-express the 
estimate (B5L) in the form (p4|), using a smaller constant k which depends only on kq and to. I 



The following estimate will be the basis of a Lyapunov function property to be established later. 
Lemma 4.2 There exist two continuous functions 

such that, for every pair of points x, z in R™ : 

(pOr) - p(z), /(*)) < -c(z)5(x, z) + (v(p(x) - p(z)), f(z)) . (26) 

Proof. As B has full column rank, there is an to x n matrix B# (for instance, its pseudo-inverse) such 
that B#B = I. We let 

v{a u ...,a n ) := ((e< fc - CTl \ . . . , e< b - CT "^)S # )' 

and 

c (C) := min e ^'^0> . 

j — l,...,m 

Now take any pair of positive vectors x, z. Denote, for each j = 1, . . . , to: 

Qj '■= (bj,p{x) - p{z)) 

and observe that 

(bi,v(p(x) - p{z))) = e qi , i = l,...,m 

so, using formula (|7|), 

m m 

(«(p(x) -p(z))J(z)) = J2J2 a ^ e{b ^ {Z)) (^'-^) ■ (27) 

i=l j=l 

Therefore (writing g{x,z) = (v(p(x) — p(z)),f{z)) for simplicity): 

m m 

(Rx)-ftz),f(x)) = ^^fl.e^^'fe-*) 

m m 

= EE a ^ e<fc3 ' Plz)>e9j (*-*) 

m m 

= ££^e< 6 -**»(e*( 9i - Qj) (e q ' - e*)) + 5 (x,z) (28) 

m m 



i=i j=i 

^ m m 

< -2C {z)Y^2a lj {q l -q ] ) 2 + g{x,z) 



i=i j=i 

~^ c o(z) Q(qi, q m ) + g(x, z) 
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where Q is the quadratic form in Lemma 4.1. Equality ( J28| ) follows by adding and subtracting g(x, z) 
and using (|27|). To justify (p9|), we note first that, for each a > 0, the function R>o — > R: 

/a(r) := e a (r - a) - e r + e a + i(r - a) 2 

is always < (because / a (0) = — e a a — 1 + e a + (l/2)a 2 < 0, f a (r) — > — oo as r — ► +oo, and /„(r) 
e a — e r + (r — a) 7^ for all r > 0). Now we use the inequality e a (r — a) — e r + e a < — \{r — a) 2 in each 
term of the sum with a = qj and r = qt (recall dy-e^'^ 2 ^ > 0). 



Lemma 4.1 gives that Q{qi, ■ ■ ■ , q m ) > k5(x, z). Thus, we may take c(z) := kcq(z)/2. 



4.1 An Entropy Distance 

Recall that we are assuming that Jq \p{r)\ dr < oo. For any fixed constant c E R, we consider the 
following function: 

Rc{r) := J p(s) ds ~ cr . 

This function is a well-defined continuous mapping R>o — > R, continuously differentiable for r > 0. 
Moreover, i? c is strictly convex, since for r > its derivative p(r) — c is strictly increasing and onto R; it 
achieves a global minimum at the unique r c £ R + where p(r c ) = c, decreases for r E [0, r c ], and increases 
to +oo for r > r c . 

The following function will play a central role: 

n 

W : R| xR^ : (x, z) » ]T R p[zi) ( Xi ) . 

i=i 

The above-mentioned properties of the functions R p ( Zi ) imply that 

x^z =>■ W(x,z)>W(z,z), (30) 

i.e., for each fixed z E R™ , the function W(-, z) has a unique global minimum, at z. Note also that the 
gradient of W(-, z): 

dW 

^(z,z) = (p(x)-p(z))' (31) 
(defined for x E R") vanishes only at x — z and that (since R p r Zj \{xi) — > +cx) if — ► +oo) 

|x| -> +oo =>■ Ly( XjZ )^ +0O (32) 
for every given z. As W(«, z) is continuous, this implies that 

{x\W(x,z) < w} (33) 

is compact for every z and every w E R. 

Remark 4.3 In the special case p = In, W(x, z) — ^ILi m 

i X i X ^ lnzj. Then this formula, 
when states x are interpreted probabilistically in applications such as chemical networks, is suggested 
by "relative entropy" considerations. □ 

4.2 Main Stability Results 

For any £ E R", let us denote by the positive class (£ + D) H^>o containing £, x^ the unique interior 
equilibrium 5 E S^, -Eq the set (possibly empty) of boundary equilibria in S^, i.e., the set S^f]Eo, and 



the set of all equilibria in S^. The main technical fact, which will imply the completeness Lemma 2.2 as 
well as the convergence parts in Theorems |l| and 0, is as follows. 
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Proposition 4.4 For each £ G R> 7 the solution x(-) of the initial value problem x — f(x) with x(Q) = £ 
is defined for alii > 0, and 

x(t) — > i? 5 as f — > +00 . 

Proo/. Fix £ in R™ . We will use 

V{x) := W(x,x^) -W(x^x^) 

as a Lyapunov-like function. By (|30|), this function is positive definite relative to the equilibrium x^, 
i.e., V(x) > for all x G R>o: an( i V( x ) = if and only if x = S*. Moreover, V is proper, meaning 
that the sublevel sets {x \ V(x) < w} are compact, for all w G R>o, by (|33|). Finally, V is continuously 
differentiable in the interior R™ , and 

VV{x)f{x) < -cS(x,x 6 ) (34) 



for all x G R™, where c = c(x^) > 0, by (31) and (Eft). In particular, since S(x,x^) = implies that x is 



an equilibrium, and since there is a unique interior equilibrium in each class: 

x G S £ n K +> & ^ ^ W(x)/(x) < . (35) 

Now consider the maximal solution x(-) starting from x(0) — £, which is a priori defined on some 
interval [0,t*). If £ G E, there is nothing to prove, so we assume that £ is not an equilibrium. As, 



in particular, £ ^ Eq, Corollary 3.13 insures that x(t) G R™ for all t > in the domain of definition. 
So y(x(t)) is differentiable for t > 0, and dV(x(t))/dt < (by @), which means that V(x(t)) is 
nondccreasing. Since V is proper, this means that the maximal trajectory is precompact, and hence it 
is defined on the entire interval [0, +oo). Furthermore, the LaSalle Invariance Theorem implies that 

x(t) — > M as ^ +oo , 

where M is the largest invariant set included in {p \ V(p) = a}, for some a > 0. 

Since x(t) G for all f, and is closed, any limit point of x(-) is included in as well. So M C 5^. 
Thus, we need only to show that M C E^. Pick any £ G M, and consider the forward trajectory z(t) 
starting from £. 

Suppose first that £ G R". Assume that ( ^ ^. Then @ says that dV(z(t))/dt < at t = 0, 
which implies that V(z(f)) < a for t small, contradicting the fact that M is included in {p | V(p) = a}. 
Thus £ = £ ■ 



Suppose now £ G R \ If C ^ ^b, then Corollary 3.13 gives that z(t Q ) G R™ for some t > 0. 
But invariance of M says that z(to) G M, and we already showed that Mf|l" C J5. So, z(to) is an 
equilibrium, which means that z(t) = z (to) and hence £ = 2 (to) ^ Rq\ a contradiction again. I 



Remark 4.5 Observe that the proof actually shows that the domain of attraction of x^ contains 
{x | V(x) < w }, where w := mm xeEo ^ ss V(x). □ 

Theorem [j] follows from the above. To prove Theorem |[ part b, pick any positive class S and £ G S. 
In the above proof, we simply note that (|35| ) holds for all x in some neighborhood of x^ = xs, so V is 
a (local) Lyapunov function for the system restricted to the class S. Finally, to see Theorem part 
c, note that if £ G Eq f] ^ then £ (being an equilibrium) is not attracted to x*, and if instead 
Eo H = then the Proposition tells us that all trajectories converge to E^ = {x^}. 



5 The Kinetic Proofreading Example, and Chemical Networks 

We now show the global stability and unique equilibrium properties of the kinetic proofreading model 
which was described in the introduction. For this purpose, we wish to see the equations as those of an 



14 



appropriate system ([j]), when restricted to a suitable class (which is determined by the constants M* 
and T*). Recalling the conservation laws M + ^2 @i = M* and T + ^ d — T*, we write the equations 
as a system of dimension n — N + 3: 

N 

4=0 

N 

M = —kiTM + k-i,iCi 

i=0 

Co = kiTM - {k-i,o + kp, ) C 



C N — kp N —\C N —\ k—i N C N . 

This is indeed a system of form (||). To see this, we use x — (T, M, Co, ■ ■ ■ , C N )' as a state, and take 
m = n- l = N + 2, 
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= 2, . . . ,to), 







ay = 0. 

Thus, D = {x | T + Co + . . . + C N = M + Co + • • • + C N = 0}, and the positive classes are of the form 
S = 5 Q /3 , intersections with M™ of the affine planes 

T + Co + . . . + C N = a , M + C + . . . + C N = (3 , 

with a > and /? > 0. The original system is nothing else than the class determined by a = T* and 
(3 = M*. Thus, the conclusions will follow from Theo rem || as soon as we prove that Sf] Eq = for any 
positive class S. To see this, we may use Proposition 3/7. Pick any x £ Mq an d an y positive class S a '@ . 
We must find some j 6 {1, . . . ,m} with the property that Xk ^ for all k G 8^ . In our application, 



Si = {1, 2} and Sj = {j + 1} for j = 2, 



If the property is not satisfied for some j £ {2, 



then Ci = for all i. But in this case, the equations for S a,t3 give that T — a > and also M = (3 > 0, 
so neither can j — - 1 be used. In conclusion, x ^ Eq, and hence part c of the theorem applies. 



5.1 Comments on Chemical Networks 

Let us very briefly indicate the motivation for the form of the systems (0) as arising from mass-action 
kinetics in chemical network theory. One studies n chemical species Ai, . . . ,A n , whose concentrations 
as a function of time are given by %i(t), . . . , x n (t) and then derives a system of differential equations for 
the Xi 's on the basis of the known reactions that occur among the substances A, . For instance, suppose 
that each molecule of A\ can react with four molecules of A 2 to produce two molecules of ^3. This is 
indicated graphically by 

A x + 4A 2 -> 2A 3 . 

Assuming that the reactor is well-mixed, the probability of such a reaction occurring, at an instant t, is 
proportional to the product x\ (t)x2(t) 4 . Since we gain two molecules of A 3 for each such reaction, this 
gives rise to a rate of increase 

±3 = 2kx\x\ (36) 
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for the concentration of A3, where k is a suitable constant of proportionality. This constant k is often 
thought of as a "reaction rate" and one writes graphically: 

A x + AA 2 ^ 2A 3 . 

In this manner, one puts together the whole system of differential equations. A convenient way to specify 
the resulting system is by building a matrix A from the reaction rates, and introducing vectors 61, . . . , b m 
to describe each of the "complexes" such as A\ + 4A 2 and 2A 3 by specifying the contributions from each 
type of molecule. For example, A\ + 4A 2 might give rise to the vector b\ — (1, 4, 0, 0, ... , 0)' and 2A 3 
to the vector b 2 = (0, 0,2,0,..., 0)', and reaction ( |36| ) contributes then the term kxix 2 (b 2 — b\) to the 
differential equations ([!]), where the reason for ll kxix 2 " is obvious, and b 2 — b\ corresponds to the fact 
that the reaction goes from the complex associated to the vector b\ to that one associated to b 2 . (Note 
how the equation for ±3 will then have the contribution kxix 2 (b 32 — b 3 i) — 2kx\x\, and also x\ will 
have a term kx\x\{b\ 2 — byy) = —kx\x\ to indicate that A\ is being eliminated at the given rate.) In 
this context, the space D is called the stoichiometric subspace associated to the reaction, and a class is 
a stoichiometric compatibility class. 

The results explained in Theorems |l| and |^ are basically the main theorems for what are called mass- 
action networks of zero deficiency and a single linkage class. We do not define these terms here. For 
more details, see for instance 0,0]. (A small remark for readers who compare our results with those 
in the chemical network literature: Condition (||) might appear to be slightly stronger than needed, 
since the zero-deficiency theorem would only require the columns bi to be affinely, rather than linearly, 
independent. However, no generality is lost, because if we start with an affinely independent set of 
vectors Vi, we can introduce the vectors bi = (1, %>[)'( i.e, just add a constant coordinate = 1) in one 
more dimension. The span of the differences m — Vj has the same dimension as D. So we need only 
consider a new set of differential equations in which we add a variable satisfying xq = 0, to bring 
everything into the setup considered in this note.) 
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